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We present a self-consistent theory for the thermodynamics of the BCS-BEC crossover in the 
normal and superfluid phase which is both conserving and gapless. It is based on the variational 
many-body formalism developed by Luttinger and Ward and by DeDominicis and Martin. Trun- 
cating the exact functional for the entropy to that obtained within a ladder approximation, the 
resulting self-consistent integral equations for the normal and anomalous Green functions are solved 
numerically for arbitrary coupling. The critical temperature, the equation of state and the entropy 
are determined as a function of the dimensionless parameter 1/kFa, which controls the crossover 
from the BCS-regime of extended pairs to the BEC-regime of tightly bound molecules. The tightly 
bound pairs turn out to be described by a Popov-type approximation for a dilute, repulsive Bose 
gas. Even though our approximation does not capture the critical behaviour near the continuous 
superfluid transition, our results provide a consistent picture for the complete crossover thermody- 
namics which compare well with recent numerical and field-theoretic approaches at the unitarity 
point. 
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I. INTRODUCTION 



The problem of a two component attractive Fermi gas 
near a resonance of the s-wave scattering length describ- 
ing the effective interaction is one of the basic many-body 
problems which has been brought into focus by the re- 
cent realization of molecular condensates in ultra-cold 
Fermi gases 0, d, ||| and the subsequent exploration of 
the crossover from a Bose-Einstein-Condensate (BEC) 
to a BCS-like state of weakly bound fermion pairs 
Clear signatures for the existence of paired fermion su- 
perfluidity with cold atoms have been provided by spec- 
troscopic measurements of the gap Q and the observa- 
tion of a vortex lattice on the BCS-side of the transition 
Q. The ability of tuning the interaction in cold Fermi 
gases through Feshbach resonances relies on the resonant 
coupling of the scattering state near zero energy of two 
colliding atoms with a bound state in a closed channel 
Q. A particularly challenging problem arises right at 
the Feshbach resonance, where the two-particle scatter- 
ing length is infinite [1, Q . Precisely at this point and 
for broad Feshbach resonances, where the range r* of the 
effective interaction is much smaller than the mean in- 
terparticle spacing [l(| [HI, [H, 0, the full many-body 
problem has the Fermi energy ep as the only energy 
scale. As pointed out by Ho \XM, the thermodynamics 
of the unitary Fermi gas is then a function only of the 
dimensionless temperature = T/Tp. More generally, 
as emphasized recently by Nikolic and Sachdev [l5|, the 
universality also extends to the behaviour away from the 
Feshbach resonance, as long as the broad resonance con- 
dition kpr* <C 1 is obeyed. Thus, for instance, the criti- 
cal temperature T c /Tp for the transition to superfluidity 
is a universal function of the inverse coupling constant 
\jkpa. 



A quantitative theoretical understanding of the many- 
body problem near a Feshbach resonance has been devel- 
oped recently through numerical calculations. In partic- 
ular, at zero temperature and for a homogeneous system, 
fixed-node Green function Monte Carlo calculations pro- 
vide quantitative results for the gap parameter (l6| . the 
equation of state [13] , and also the momentum distribu- 
tion, the condensate fraction and the pair size [IH of the 
ground-state for arbitrary values of l/kpa. As expected 
in the case of an s-wave resonance [19J, these quanti- 
ties all evolve continously as the coupling is varied from 
the BCS to the BEC-limit. An important ingredient in 
these results is their account for the repulsive interaction 
between strongly bound dimcrs in the BEC-limit with 
scattering length add ~ 0.60 a > [2fi| , This interac- 
tion is missing in the early qualitative descriptions of the 
T = BCS-BEC crossover problem by Eagles [2l| and 
Leggett [2a |, which are based on using the standard BCS- 
groundstate as a variational Ansatz for arbitrary coupling 
|23j. Beyond a purely numerical approach, the BCS- 
BEC crossover problem has recently become amenable 
also to analytical methods via an e = 4 — d expansion 
p3 |. It is based on the observation [25| that at the uni- 
tarity point in d = 4 (i.e. the point where a two-particle 
bound state appears) the two-component Fermi gas is in 
fact an ideal Bose gas, because a zero range interaction 
in d = 4 can bind a state only at infinitely strong at- 
traction. In two dimensions, in turn, binding appears 
at arbitrary small couplings and the unitar y F ermi gas in 
d < 2 coincides with a non- interacting one [25| . Within a 
field theoretic description, the physically interesting 3D 
problem can thus be approached by extrapolating ex- 
pansions from the upper and lower critical dimensions 
d = 4 and d = 2 respectively [26]. At finite temper- 
ature, numerical calculations arc available for the ther- 
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modynamics at the unitarity point. They are based on 
an auxiliary field quantum Monte-Carlo method for the 
continuum problem [27| and on a diagrammatic determi- 
nant Monte-Carlo method for the negative [/-Hubbard 
model [28| . Field-theoretic results at finite temperature, 
which open the possibility for controlled and systematic 
expansions for the crossover thermodynamics have been 
obtained very recently by Nishida [26| within an expan- 
sion around both the upper or lower critical dimension 
and by Nikolic and Sachdev [l5| within &1/N expansion 
for a 2N component Fermi gas. 

Our aim in the following is to present a self-consistent 
many-body theory for the thermodynamics of resonantly 
interacting fermions at arbitrary temperatures and de- 
tuning, which directly addresses the physically relevant 
case of a three dimensional, two-component Fermi gas. 
The theory is based on a conserving, so-called ^-derivable 
approach to the many-body problem, in which the exact 
one- or two-particle Green functions serve as an infinite 
set of variational parameters. It is an extension of earlier 
work by one of us 0, [3(| HH and employs a combina- 
tion of the Luttinger-Ward and DeDominicis-Martin ap- 
proach for obtaining the grand canonical potential and 
the entropy, respectively. The condition of gaplessness is 
enforced by a modified coupling constant, thus account- 
ing for the proper low energy behaviour in terms of a 
Bogoliubov- Anderson mode. We provide quantitative re- 
sults for the critical temperature, the equation of state 
and the entropy near the Feshbach-resonancc as a func- 
tion of both T/Tp and 1/kpa. In spite of the fact that the 
critical behaviour at the continuous superfluid transition 
is not captured correctly in our approach, which gives 
rise to a weak first order transition, the results provide a 
quantitative and consistent picture of the crossover which 
obeys thermodynamic relations at the percent level. Our 
variational method is complementary both to purely nu- 
merical and to field theoretic approaches to the problem. 
The results can be used e.g. to predict the final tempera- 
ture reached after an adiabatic ramp across the Feshbach- 
resonance starting deeply in the BEC-regime [HJ or to 
determine the size of the atom cloud in a harmonic trap 
near unitarity as a function of temperature. 

The paper is organized as follows: in Sec. |TT] we in- 
troduce our model and the basic many-body formalism 
necessary for deriving a set of self-consistent equations for 
the Green and vertex functions which are the variational 
parameters of the theory. The complete thermodynam- 
ics is then determined by integrals of the momentum and 
frequency dependent Green functions. It is shown that 
with a modified coupling constant, the theory can be for- 
mulated in a way consistent with Ward identities, which 
guarantees a gapless Bogoliubov- Anderson mode for arbi- 
trary strength of the coupling. In Sec. IIIII we discuss the 
numerical solution, providing quantitative results for the 
critical temperature, the pressure, internal energy and 
the entropy of the BCS-BEC crossover both in the nor- 
mal and superfluid phase. They are compared both with 
experimental and theoretical results based on numerical 



and field-theoretic approaches. Finally in Sec. lIVI wc give 
a brief summary, and indicate open problems. 

II. A MANY-BODY THEORY OF 
RESONANTLY INTERACTING FERMIONS 

In order to describe interacting fermions near a Fcsh- 
bach resonance, it is in general necessary to include 
the resonant, closed channel bound state explicitely, e.g. 
within a Bose-Fermi-resonance model [33|, H3- As has 
been shown for instance by Diener and Ho [111 ], how- 
ever, the situation can be simplified in the case of broad 
Fcshbach resonances, where the effective range r* of the 
resonant interaction is much smaller than both the back- 
ground scattering length a.b g and the Fermi wavelength 
Xp. In this limit, which is in fact appropriate for the 
existing experimental studies of the BCS-BEC crossover 
problem in 6 Li [j| and in 40 K [l|, the problem can be 
reduced to a single channel Hamiltonian with an instan- 
taneous interaction [lj], HH [H, EH ■ The associated effec- 
tive two-body interaction is thus described by a pseudo- 
potential V(r) ~ 8(r) (appropriately renormalizcd, see 
below) with a strength proportional to the scattering 
length 

« = a fc3 (l-^y. (2.1) 

Here a,b g is the off-resonant background scattering length 
in the absence of the coupling to the closed channel while 
AS and B describe the width and position of the reso- 
nance which may be tuned by an external magnetic field 
B. The interacting Fermi system is thus described by the 
standard Hamiltonian 

H=[ d d r£|^[V^(r)][Wv(r)] 

(7 

+ ]-Jd d rJd d r'Y,V(r-v>) (2 ' 2) 

xi + (r)i+(r')i,(r')i(r) , 

where tp a (r) and ip£ (r) are the usual fcrmion field opera- 
tors. The formal spin index a labels two internal degrees 
of freedom, which in practice are two different hypcrfine 
states. In the approximation, where the effective range of 
the resonant interaction is taken to zero, the interaction 
potential can formally be replaced by a delta potential 
between fermions of opposite spin 

V(r-r')=g Q S(r-v'). (2.3) 

Its strength go needs to be renormalizcd for dimensions 
d > 2 by introducing the scattering amplitude g via 

1 _ 1 f d d k m 

~ 9 = 7o + J J^y m? ' {2A > 

For dimensions d > 2 the integral diverges at high mo- 
menta. Since the scattering amplitude g is kept constant, 
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the bare interaction parameter go must be taken to zero 
in the limit where the cutoff diverges. The associated 
limiting process go — ► —0 accounts for the replacement 
of the bare potential (|2.3[) by a pseudopotential with the 
proper scattering length. While the formulas are derived 
for arbitrary space dimensions d, eventually we consider 
fermions for d = 3. In this case the scattering amplitude 
g is simply connected to the s-wave scattering length a 
given in (|2.1[) by g = 4Trh 2 a/m. 

In the following, we consider a homogeneous situation 
described by a grand canonical distribution at fixed tem- 
perature and chemical potential. The thermodynamic 
properties thus follow from the grand partition function 

Z = Tr{exp(-0[H - fiN])} (2.5) 

and the associated grand potential 

il = n(T,(i) = -p- 1 ]nZ (2.6) 

which is directly related to the pressure p via Q = —pV. 
Within our simplified model, where the range of the in- 
teraction is set to zero, the Fermi system is described 
by three parameters: the temperature T, the chemical 
potential fi and the s-wave scattering length a. Apart 
from an overall scale, the thermodynamics thus depends 
only on two dimcnsionless ratios. It is convenient to re- 
place the chemical potential \i by the fermion density 
n = k F /3ir 2 , which defines the Fermi wave number Uf 
and the Fermi energy £f = h 2 k F /2m as characteristic 
length and energy scales. The equilibrium state is then 
uniquely determined by only two parameters: the dimen- 
sionless temperature 8 = T/ep (we choose units for the 
temperature in which ks = 1), and the dimcnsionless in- 
teraction strength v = 1/kpa. In the special case B = Bo 
of an infinite scattering length (the so-called unitarity 
limit), the parameter v drops out and the resulting ther- 
modynamic quantities are universal functions of 9 [Til ]. 

A. Luttinger-Ward formalism 

The BCS-BEC crossover is controlled by two physical 
phenomena. The first one is connected with the forma- 
tion of pairs due to the attractive interaction. The sec- 
ond one is the transition to superfluidity below a certain 
critical temperature T c . In the BCS-limit, the formation 
of pairs and the superfluid transition are simultaneous. 
The transition is driven by the thermal breakup of pairs, 
i.e. by excitations which may be described by a purely 
fermionic theory With increasing strength of the inter- 
action, however, there is an increasingly wide range of 
temperatures where bound pairs coexist with unpaired 
fermions. In the BEC-limit, pair formation, as a chem- 
ical equilibrium between bound and dissociated atoms, 
occurs at a temperature scale much higher than the su- 
perfluid transition. The latter is driven by collective ex- 
citations of a then purely bosonic system. A proper de- 
scription of the crossover thus requires to account for 
both bosonic and fermionic excitations simultaneously. 



Following the formalism developed by Luttinger and 
Ward [IH for non-superfluid interacting Fermi systems, 
the grand thermodynamic potential (|2.6|) can be ex- 
pressed as a unique functional of the Green function 

G^{v-v',t-t') 

_( S aa ,G(r-r',T-r') e a(J >T{v - r',r- r') \ 
\-e ac .T*{v> -r,r- r') -<W£(r' - r, r' - r) J 

(2.7) 

in the form 

ft[G] =/?- 1 (-iTr{-lnG+[G - 1 G-l]}-$[G]) . (2.8) 

The trace Tr is defined with respect to the formal index 
X = (r, t, a, a) which combines the space variable r, the 
imaginary time r, the spin index cr, and the Nambu index 
a. The interaction between the fermions is described by 
the functional $[G], which can be expressed in terms of 
a perturbation series of irreducible Fcynman-Diagrams 
where the propagator lines are dressed and identified by 
the matrix Green function G of (|2.7p . 

While the formalism of Luttinger and Ward was origi- 
nally derived for normal quantum liquids, it is well suited 
also to describe superfluid systems. Indeed the nondi- 
agonal elements of the matrix Green function G rep- 
resent the order parameter of the superfluid transition. 
The minimization of the grand potential Q[G] as a func- 
tional of the Green function G thus incorporates the stan- 
dard thermodynamic criterion that the order parame- 
ter is found by minimizing the thermodynamic potential. 
The stationarity condition 

Sn[G]/8G = (2.9) 

uniquely determines the full matrix Green function G of 
the interacting system and hence the order parameter. It 
is important to note, that the thermodynamic potential 
fl[G] depends on the exact Green function G. The for- 
malism of Luttinger and Ward thus leads via (|2.9[) to a 
self-consistent theory for the matrix Green function G. 
Since the Green functions contain information about the 
full dynamical behaviour via the imaginary time depen- 
dence of the Matsubara formalism, the Luttinger-Ward 
approach not only provides results for the equilibrium 
thermodynamic quantities but also determines spectral 
functions and transport properties. In our present work, 
however, dynamical properties will not be discussed. 

The functional $[G] is defined by an infinite perturba- 
tion series of irreducible Feynman diagrams and an exact 
expression for $[G] is clearly beyond what can be done 
analytically. An approximation which properly describes 
the formation of pairs, is a ladder approximation [36| . In 
Fig. [TJ the related diagrams of $[G] are shown. The lad- 
der approximation is self consistent because the propaga- 
tor lines are dressed lines which are identified by the ma- 
trix Green function G. In the weak coupling BCS regime 
the ladder approximation becomes exact. For very strong 
attractive interactions, well above the pairing threshold, 
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FIG. 1: The functional $[G] in self-consistent ladder approx- 
imation. The propagator lines are dressed lines identified by 
the matrix Green function G. 



the fermion system is a Bose liquid of dilute atom pairs. 
In this limit the ladder approximation describes the for- 
mation of pairs (two-particle problem) exactly, however 
the interaction between the pairs (four-particle problem) 
only approximately 0, [3?| . In particular the resulting 
dimer-dimer scattering length is given by the Born ap- 



proximation a 



(B) 



2a. 



As in standard thermodynamics, the entropy (|2.10|) 
is maximized under the constraints that all conserved 
quantities are kept constant. For the interacting fermion 
system defined by the Hamiltonian (|2.2p the conserved 
quantities are the internal energy U = (H) and the par- 
ticle number N = —^Ti{G}. Evaluating the thermal 
average of the Hamiltonian (|2.2p we find that U can be 
expressed in terms of G and T (see (|2.15j) and (|2.24|) be- 
low). 

Consequently, the entropy [G, Y] = S[G, T], the in- 
ternal energy U[G, T], and the particle number N[G] arc 
functionals depending on G and T. In order to find the 
maximum of the entropy under the constraint of given 
average values of the particle number and the internal 
energy, DeDominicis and Martin [38[ consider the func- 
tional 

w[G, r] = [G, r] - \uU[G, r] - \ n n[g] (2.11) 



B. DeDominicis-Martin formalism 

An extension of the Luttinger-Ward formalism was 
given by DeDominicis and Martin (38[. They introduce 
up to four external fields, which couple to products of 
one, two, three, and four field operators, and perform 
the Legendre transformations to the corresponding con- 
jugate variables - the Green functions. For fermion sys- 
tems only two external fields arc relevant which couple to 
even products of fermion field operators. The related two 
conjugate variables of the Legendre transformation are 
the one-particle Green function G and the two-particle 
Green function G 2 . Within our approach below, the 
second Legendre transformation is performed explicitely. 
A more convenient conjugate variable is then the ver- 
tex function V which is related to G2 by (|2.15|) below. 
Thus, DeDominicis and Martin obtain a thermodynamic 
potential which is a functional of both G and V. More 
precisely, it turns out that the relevant functional is the 
entropy S = F^ where 

FW[G,T\ =iTr{-lnG+ [(-itiu n )G - 1]} 

+ iTr{ln[l - if] + if + i[if] 2 (2.10) 
- (l/4!)[f] 2 } + /C (2) [G,r] 

(see (61) in the second paper of Ref.[H and identify G\ = 
G, G 2 = -r, and G 2 = -f therein), f is defined in (|2T3|) 
below. 

The formalism of DeDominicis and Martin is ideally 
adapted to describe the BCS-BEC crossover because it 
explicitly deals with the one-particle Green function G, 
which represent the properties of the single fermions, 
and the vertex function Y, which describes the eventu- 
ally purely bosonic properties of the fermion pairs (both 
condensed or noncondensed) . In particular, a full imple- 
mentation of their formalism is needed to correctly ac- 
count for four particle correlation, which is necessary to 
obtain the exact result add = 0.60 a for the dimer-dimer 
scattering length in the BEC limit. 



where Xu and Xn are two Lagrange parameters for the 
two constraints. Alternatively and equivalently, we con- 
sider the functional 



n[G, r] = u[G, r] — t s[g, r}-^ n[g] 



(2.12) 



which is the grand thermodynamic potential where the 
temperature T and the chemical potential \i are the La- 
grange parameters. Both functionals (|2.1ip and (|2.12[) 
must be stationary under small variations of G and T. In 
this way, we obtain the stationarity criteria 



Sn[G,T}/6G = and SQ.[G,Y]/SY = 



(2.13) 



which uniquely determine the one-particle Green func- 
tion G and the vertex function T. 

In order to simplify the second trace in the entropy 
functional (|2 . 10[) it is convenient to define a modified ver- 
tex function f by 



XiX 2 X 3 X^ 



— U X 1 Y 1 '~ r X2Y 2 l Y l Y 2 Y 3 Yi f ~'YzX 3 U ^ 



r Y 4 X 4 

(2.14) 

where the four external propagator lines are amputated 
only half way (see (46) in the second paper of l3o). For 
a proper definition of the second trace and the related 
matrix products the four indices of the modified ver- 
tex function must be grouped into pairs according to 
f = f (x 1 x 2 )(x 3 x 4 )- The last term in (|2.10[) . the func- 
tional /C( 2 )[G,r] (depicted in Fig. [5]), is defined by an 
infinite perturbation series of 2-line irreducible Feynman 



/c( 2 )[G,r] = 





FIG. 2: The functional /C (2) [G, F] is the sum of all 2-line irre- 
ducible diagrams. The propagator lines and the vertices (full 
circles) are dressed and identified with G and T, respectively. 
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diagrams, where the propagator lines and the vertices are 
dressed and identified by the one-particle Green function 
G and by the vertex function T, respectively. 

In order to understand the physical meaning of the 
various contributions to the thermodynamic potential, 
we note that the Luttinger-Ward formalism and the 
DcDominicis-Martin formalism are related to each other 
by a Legendre transformation, in which the bare two- 
particle interaction as an external field is transformed 
into the two-particle Green function G2. This Legendre 
transformation may be interpreted as a renormalization 
procedure. Since the two-particle Green function Gi is 
expressed in terms of the vertex function T by 

G2.x 1 x 2 x 3 x 4 =Gx 1 x 3 Gx 2 x 4 

— Gx 1 x 4 Gx 2 x 3 — GxxXiGxsXi 

- Gx 1 Y 1 Gx 2 Y 2 ^Y 1 Y 2 Y 3 Y i GY 3 X 3 GY i X 4 , 

(2.15) 

the bare two-particle interaction is replaced by the vertex 
function T as the renormalized interaction, which is the 
many-particle generalization of the scattering amplitude 
9- 

We may compare the functionals (|2.8p and (|2.10p di- 
rectly with each other. The first trace in (|2. 10[) is iden- 
tified by the trace in (|2.8[) . which describe the contribu- 
tion of single particles to the grand canonical potential 
and entropy, respectively. The second trace and the func- 
tional JC^ [G, r] in (|2.10p which represent the interaction 
terms are related to the functional <f> [G] in (|2.8p . By close 
inspection (see (62) in Ref. [H) we find that the second 
trace in (|2.10p represents the inverted perturbation series 



of ladder diagrams. It includes both particle-particle and 
particle-hole ladders, which describe the scattering and 
formation of pairs and it also includes bubble diagrams, 
which describe the screening of the interaction. 

These three types of diagrams and mixtures of them 
arise because the vertices are symmetrized so that each of 
them can be expressed as a sum of three unsymmetrized 
vertices. As a result, the self-consistent ladder approx- 
imation of the functional $[G] shown in Fig. [1] can be 
reformulated within the formalism of DcDominicis and 
Martin in the following way: the second trace is ap- 
proximated by keeping only the particle-particle ladders 
and the complicated functional /C' 2 '[G, T] is set equal 
to zero. This neglects the screening of the interaction 
due to particle-hole excitations (see Subsec. IIII Al below) 
and also the coupling between collective excitations and 
bound pairs. 

In the following subsections we employ the formal- 
ism of DeDominicis and Martin to construct explicit ex- 
pressions for S[G,T], U[G,T], and N[G}. From (l2~12l 
we obtain the functional fi[G,F]. The stationarity cri- 
teria (|2.13p imply two self-consistent equations for the 
Green function G and the vertex function T. Solv- 
ing the second equation with respect to T and insert- 
ing the resulting vertex function into S1[G, T] wc recover 
the functional fl[G] of the Luttinger-Ward formalism to- 
gether with the stationarity condition (|2.9p . This fact 
explicitely demonstrates the equivalence of the Luttinger- 
Ward and DcDominics-Martin formalism for our approx- 
imation scheme (see (|2.41l) - (|2.43l) below) once the appro- 
priate stationarity conditions have been taken into ac- 
count. 



C. Thermodynamic potentials 

The formalism of Luttinger and Ward [35[ allows to calculate directly the grand thermodynamic potential O. The 
functional ^[G] has been evaluated explicitly in Ref. [3l|. Inserting this result into (|2.8p we obtain 

d d k 1 



n{G]=-L d J I^Tr{-ln[G(k,^)] + [G (k,c n )- 1 G(k, W „)-l]} 

+ L d g,\T{0M 2 + \L d J 0^- d i^Tr{ln[l+.g oX (K,aO]} • 



In this formula the matrix Green function is defined by 



Knowledge of the matrix Green functions determines the matrix pair propagator via 

/ f d d k 1 v-^ 

X(K, fi„) = (Xaa> (K, fi„)) = J^fjd p Z_] Gaa ' ( K ~~ k ' fi ™ ~ u n) G °"*' ( k ) ^n) 



(2.16) 



G(k,.„) = (G aa ^ n) ) = • (2.17) 



(2.18) 



In order to distinguish between fermionic and bosonic functions the fermionic wave vectors and Matsubara fre- 
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quencies are denoted by small letters, while the bosonic 
wave vectors and Matsubara frequencies are denoted by 
capital letters. In the second term of (|2.16j) the anoma- 
lous Green function is identified by JT(0,0) = J-(r = 
0, r = 0). The formulas are derived for an arbitrary di- 
mension of space d. The volume is assumed to be a cube 
with edge length L and periodic boundary conditions, 
where the limit L — > oo is taken. 

The strength of the attractive interaction is included 
by the bare interaction parameter g . The kinetic energy 
of the atoms = fi, 2 k 2 /2m and the chemical potential 
fi are implicitly included via the free matrix Green func- 
tion Go(k, u n ), which is related to the free normal Green 
function 

S (k,O = l/[-ifej n + ek - (2.19) 
and the free anomalous Green function 

.Fo(k,w„) = (2.20) 
by a formula which is analogous to (|2.17[) . The tem- 



perature T is included explicitly by the factors 1//3 and 
implicitly by the Matsubara frequencies u„ and fl n . 

As evident from (|2.16|) . the formalism of Luttinger and 
Ward, though including the exact single-particle Green 
function, still contains the bare coupling constant go. In 
the DcDominicis-Martin formalism the bare coupling is 
renormalized and replaced by the exact vertex function T 
via a second Legendre transformation. The correspond- 
ing functional S[G, T] = F 2 [G,T], is just the dimension- 
less entropy as given in f (|2.10p ) 

As discussed above we restrict the second trace in 
(|2.10j) to the particle-particle ladders and by the nature 
of our interaction potential (|2.3|) to s-wavc scattering. 
Furthermore, we omit the 2-line irreducible Feynman dia- 
grams by setting K,^ [G, T] = 0. This approximation cov- 
ers the essential features of the crossover problem namely 
the formation of pairs and their condensation. Within 
our ladder approximation, the DcDominicis-Martin for- 
malism thus leads to an expression for the entropy of the 
form 



S[G,T] =PL d J -0^ i^Tr{-ln[G(k, W „)] + H^„G(k,c„)-l]} 

+ \pL d J 0^ i^Tr{ln[l-x(K,fi„)r(K,O n )]+x(K,f2„)r(K,0„)} . 



(2.21) 



The first term is clear. It is directly obtained from the 
first trace in (|2.10p . However, the second term resulting 
from the second trace in (|2.10p needs further explanation. 
From (|2.14p and the definition of the pair propagator 
(gUBH we infer 

f (k, n n ) = x (k, n„) 1/2 r(K, n n ) x (k, a,) 1 / 2 

= x(K,n n )r(K,n n ) (2.22) 

= r(K,Q„)x(K,0„) . 

The reduction to particle-particle ladders implies that 
the Nambu indices are pairwise identical. In this way, 
the four Nambu indices of the vertex function T reduce 
to two Nambu indices. As a result, the vertex function 
T(K, fl n ) = (r QQ /(K, f2 n )) is a 2x 2 matrix in the Nambu 
space similar to the matrix Green function (|2.17p . For 
the formalism of Luttinger and Ward the reduction of the 
vertex is described in detail in Ref. [U and also in Ref. 
|29| . Since the second trace of (|2.10p is reduced to the 
particle-particle ladders and the structure of the vertex 
function is simplified considerably due to s-wave scatter- 
ing, the prefactors of the terms in the second trace of 
(|2.2ip are changed. The factor ^ in front of f disap- 
pears. Furthermore, the quadratic terms in the second 
trace cancel. 



Another important quantity to consider is the internal 
energy U. With the help of the delta potential (|2.3p we 
find for the expectation value of the Hamiltonian (|2.2p . 

J 2m 

a 

+ IJ ^E-9o(^( r )V^W^(r)^(r)) . 

GO 1 

(2.23) 

The second term contains an average of four fermion field 
operators which can be expressed in terms of the two- 
particle Green function G2. Following the formalism of 
DeDominicis and Martin [38| and using (|2 . 15[) the two- 
particle Green function can be expressed by four terms. 
The first three terms represent the three possibilities to 
factorize the two-particle Green function into products of 
two one-particle Green functions according to the Wick 
theorem. These terms provide the Hartree energy, the 
Fock energy, and the Bogoliubov energy. The fourth term 
is the connected part of the two-particle Green function 
and provides the correlation energy. Taking all terms 
together we obtain the internal energy 
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U[G,T] 



2L d 
1 



d d k 



1 s k g(k,T = -o) + L d 9a \T{o,Q)Y 



(27T) 

d d JC 1 



L d 

2 ./ (27r) d (3 



J2 9o Tr{ X (K, O n ) - X (K, ft„) r(K, ft n ) X (K, 0„)} 



(2.24) 



Finally, the particle number N is defined by the aver- 
age 



N = (N) 



d d r 



5>J(r)V») 



(2.25) 



which can be expressed in terms of the normal Green 
function in the standard form 



N\G] 



-2L C 



d d k 



(27T) 



3 S(k,r = -o) 



(2.26) 



The entropy (|2.21|) . the internal energy (|2.24p . and the 
particle number (|2.26| are the basic functionals of the for- 
malism of DcDominicis and Martin. The entropy S[G, T] 
is maximized under the constraints that the internal en- 
ergy U[G, r] and the particle number N[G] are constant. 
In order to do this, the grand thermodynamic potential 
fl[G, F] is defined by (|2.12p where the temperature T and 
the chemical potential /j, are Lagrange parameters. The 
self-consistent equations for the Green function G and 
the vertex function T are obtained from the stationar- 
ity conditions (|2.13[) . Formally, the formalism of DcDo- 
minicis and Martin yields a different expression for the 
grand thermodynamic potential f2 than the formalism of 
Luttinger and Ward does by (|2.16|) . However, it can be 
shown that the results are identical if and only if G and 
r satisfy the self-consistent equations (see end of Subsec. 

MM- 

The functionals (|2~2"Tj) , (|2~24"p . and (|2~2cj|) do not de- 
pend explicitly on the thermodynamic parameters T and 
(i. While the temperature appears explicitely via the fac- 
tor (3 = 1/T and the Matsubara frequency u) n ~ Q n ~ T, 
a proper rescaling of the functions G^>(3G, x^PX: 
and r — > r implies that all factors (3 and T cancel in all 
three functionals. The temperature T and the chemical 
potential /x enter only as Lagrange parameters via the 
constraints. This fact is a general property of the for- 
malism of DeDominicis and Martin. The fermion mass 
m, the kinetic energy £k = fi, 2 k 2 /2m, and the interaction 
parameter go , which determine the microscopic proper- 
ties of the interacting fermion system, are present only 
in the internal energy functional (|2.24|) . 

An alternative expression for the entropy is obtained 
from the grand thermodynamic potential of Luttinger 
and Ward (|2.16p according to the standard thermody- 
namic relation 



Taking the partial derivative we obtain an expressions 
which formally differs from (|2.2ip . However, provided 
that G and V satisfy the self-consistent equations, the 
results for the entropy will be identical. Therefore, both 
the Luttinger- Ward and the DcDominicis-Martin formal- 
ism exactly obey all the standard thermodynamic re- 
lations provided the Green functions obey the station- 
arity conditions (|2.9[) and (|2.13[) . The equivalence of 
the different formal expressions in thermal equilibrium is 
very important for the consistency of our theory and the 
compatibility of the self-consistent ladder approximation 
for all thermodynamic quantities. Apart from the en- 
tropy, we can also determine the pressure p — —£l/L d 
as a functional of the Green function G using (|2.16p 
or (|2.12p . The dimensionless thermodynamic quantities 
£}/Nsf, U/Nef and S/N will be calculated numerically 
in Sec. IIIII and discussed in the following sections. 



D. Self-consistent equations for the Green and 
vertex functions 



The self-consistent equations for the Green functions 
follow directly from the stationarity condition (|2.9[) . In- 
serting the general functional of the Luttinger- Ward for- 
malism (|2.8|) into this condition we obtain the Dyson 
equation 

G- a l (k, u n ) = GoX* (k, a/») - (k, u n ) . (2.28) 
The self energy E is identified by the functional derivative 



1 ggjg] 



(3L d 5G a , a (k,u n ) 



(2.29) 



The functional 3>[G] is defined by a perturbation series. 
The related Feynman diagrams are shown in Fig.[T]for the 
self-consistent ladder approximation. Inserting the grand 
thermodynamic potential (|2.16| into the constraint (|2.9|) . 
we obtain an explicit expression for the self energy which 
is 



E aa '(r,r) =Ei, QQ ' 8(y)5 f {t/K) 

+ G a > a {-r, -r)r Qa /(r,r) 



(2.30) 



In the first term 8f{t/K) is the fcrmionic delta function 
which is antipcriodic. The order parameter of the super- 
fluid transition 



S = -dVL/dT 



(2.27) 



A = . go ^(0,0) 



(2.31) 
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is represented by the nondiagonal elements of the matrix 

I A* o) ■ (2 ' 32) 



In the second term of (|2.30p T is the matrix vertex func- 
tion, which is related to the matrix pair propagator \ 
by 



T Q 1, (K, fi„) = g Q 1 S aa > + Xaa> (K, tt n ) 



(2.33) 



Eventually, x ' IS represented in terms of the matrix Green 
function by (f2TT8|) . Eq. (|2~33)) is just the Bethe-Salpeter 
equation in ladder approximation. It is responsible for 
the fact that the binding of fermion pairs is described 
appropriately. Taken together, we have now a set of self- 
consistent equations for the matrix Green function G and 
the matrix vertex function T which have to be solved 
numerically. 



Alternatively we can derive the self-consistent equa- 
tions by inserting the functional (|2.12|) of the formalism 
of DcDominicis and Martin into the related stationarity 
conditions (|2.13[) . We obtain the Dyson equation (|2.28p 
from the first condition and the Bethe-Salpeter equa- 
tion (|2.33|) from the second condition. In this way we 
prove that both the Luttinger-Ward formalism and the 
DcDominicis-Martin formalism are equivalent within our 
approximation. 



Unfortunately, in the present form, the matrix pair 
propagator x defined in (|2.18|) is divergent. While the 
sum over the Matsubara frequencies is finite, the integral 
over the wave vector is ultraviolet divergent for dimen- 
sions d > 2. For this reason a renormalization is neces- 
sary. We define the regularized pair propagator by 



M aa ,(K,n n ) = J ^A[I^G aa ,(K-k,O x 



u n )G aa i (k, uj n ) 



(2.34) 



r 



Inserting this formula into (|2.33[) we obtain the renor- 
malized Bethe-Salpeter equation 



I 1 -*, (K, Q n ) = g- l 8 aa , + M aa , (K, il r 



(2.35) 



The bare interaction strength go is renormalized accord- 
ing to (|2.4p and replaced by the scattering amplitude g. 
For d = 3 dimensions g is expressed in terms of the s-wave 
scattering length a by g = 4irh 2 a/m. 

The zero range of the interaction between the fermions 
implies that 



(2.36) 



is infinite. For this reason the order-parameter formula 
(|2.3ip must be renormalized, too. Replacing the bare 
interaction strength go by the scattering amplitude g ac- 
cording to (|2.4p , we obtain the renormalized formula 



A 



d d k r 



(27r) d 



.F(k, t = 0) + A 



m 



(2.37) 



Here, the integral over the wave vector is finite. 



E. Reformulation in terms of mean-field Green 
functions 

In mean-field approximation the self energy S(k, uj n ) is 
replaced by £i defined in (|2.32p . Since £i depends nei- 
ther on wavevector nor on frequency, the approximation 



£ f=a Ei just describes the formation of a pair condensate 
within a BCS-type mean-field theory where the destruc- 
tion of superfluidity is driven by the breakup of pairs. 
This is the correct description in the weak coupling limit, 
however for strong coupling the superfluid transition is 
driven by finite momentum pairs whose contribution is 
contained in the second term of the self energy (|2.30p . 
Inserting the mean-field self energy into the Dyson equa- 
tion (|2.28p we obtain 



G^k,^)" 1 =G (k, w „)- 1 -£ 1 

-iu) n + (e k - ^) -A 

-iuj n ~ (e k - fx) 



-A* 



(2.38) 



where Gi is the matrix Green function in mean-field ap- 
proximation. 

If we consider the self-consistent equations and the for- 
mulas for the thermodynamic potentials we realize that 
the spectrum e k of the fermionic atoms and the chemi- 
cal potential [i enter the formulas only implicitly via the 
free matrix Green function Go- We can transform the 
formulas so that Go is replaced in favor of the mean-field 
matrix Green function Gi. As a result we obtain the 
Dyson equation 



G _1 ,(k,w„) 



-(k,u; n ) (2.39) 



where 



Eaa' ( r , t) — G a ' a (— r, — r)T aa i (r, r) 



(2.40) 
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is the second term of the self energy (|2.30[) . The other thermodynamic potential (]2 . 16[) is transformed into 
self-consistent equations remain unchanged. The grand 



Q = ~ Ld J Wf ^^{-HG{^n)] + [G^ n )-^G{h,u> n )-l]} 

\ Ld j 0f\T.n-Hn^n)/ g ,]} . 



- L d — 



(2.41) 



50 

For a combination of the internal energy ()2.24j) and the particle number (|2.26| we obtain the formula 
U-[iN = -L d J -0j i^Ti-{[G 1 (k,u..„)" , +irw,]G(k,a,„)} 

- i ^-5 ii /if|^ T ' {|r(K ' !i - )/s "- 11} - 

The entropy (|2.21|) depends neither on Go nor on G\ . We transform the formula into 
5 =P L<1 J Wj* \ ^ Tl "^ ln[G(k ' Un)] + H^«G(k, u n ) - 1] } 

-\l* Ld J 0f ^E^{-Hr(K,o„)/ 5 o] + [r(K,f2„)/ ffo -i]} . 



(2.42) 



(2.43) 



In the above three formulas we have simplified the terms involving the vertex function T by using the Bethe-Salpeter 
equation (|2.33[) . The grand thermodynamic potential (|2.41[) was derived using the formalism of Luttinger and Ward 
35f while the other two quantities (|2.42p and (|2.43p were derived using the formalism of DcDominicis and Martin 
381 ]. It is now not hard to see that the above expressions obey the thermodynamic relation 

n = U -TS - fj,N (2.44) 

which cxplicitely shows that both formalisms arc indeed equivalent yielding the same results for all thermodynamic 
potentials in self-consistent ladder approximation provided G and Y satisfy the appropriate stationarity equations. 

F. Mean- field approximation 

If we insert G\ for the matrix Green function G and neglect all terms containing the vertex function T we obtain the 
thermodynamic potentials in mean-field approximation. In particular the mean-field grand thermodynamic potential 
is given by 



Or = -L*j i£Tr{-ln[G l( k jW „)]} -L*. 



9a 



(2.45) 



while the mean-field formula for the combination of the internal energy and the particle number arc 

U 1 - t xN 1 = -L d / 44 ^Y t Tr{[G 1 Oc,u n )- 1 + itiu n ]G 1 (k,w n )}-L d ^- , (2.46) 
J (27r) d P j£ 9o 

and the mean-field entropy is 

Si=(3L d J -0^ i Tr i - ln t Gl ( k ' ^ + ( k > w») " 1]} ■ ( 2 -47) 
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In order to obtain finite results, we must define the sums over the Matsubara frequencies as described in Appendix 
IA"1 The sums can be evaluated explicitly. This yields 

ill = E -j 3 2L d J l n [i + cxp(-/3(£ k - n)] , (2.48) 

Ui - pNi = E Q + 2L d J (E k -n) n k , (2.49) 

Si = -2L d J -^A{(l-n k )ln(l-nk) + n k hin k } , (2.50) 

which are the well known results of a BCS variational Ansatz for arbitrary coupling where 

\Al 
9o 

I 



E o = -2L d J-0^ i[(£ k - M )-(e k - M )] - L d - 



(2.51) 



is an energy constant which after renormalization go — > 
g (see (|2.57p ) reduces to the BCS condensation energy. 
Here -E k is the spectrum of the quasiparticles, defined by 

(i? k -M) = [( £k -/i) 2 + |A| 2 ] 1 / 2 , (2.52) 

and n k denotes the Fermi distribution function of the 
quasiparticles 

7i k = l/[exp(/3(£ k - n)) + 1] . (2.53) 

We find that the regularization of the Matsubara-fre- 
quency sums described in Appendix [A] affects only the 
energy constant E . The regularization has been cho- 
sen such that for zero interaction the results for the ideal 
Fermi gas are obtained which implies Eq = 0. The other 
terms in (|2.48|) - (|2.50[) are not affected by the regulariza- 
tion. 

G. Beyond mean-field 

In the mean-field approximation, the formation and 
condensation of fermion pairs occur at the same temper- 
ature. This is the well known BCS scenario, which is 
perfectly captured by the exact solution of the reduced 
BCS-Hamiltonian. Formally, this solution can easily be 
extended to arbitrary coupling strengths [23}. At zero 
temperature, it provides a smooth crossover from the 
BCS groundstatc of highly overlapping pairs to a perfect 
Bose-Einstein condensate at infinite coupling, similar to 
the variational Ansatz of Eagles and Leggett [2l|, [22] • At 
finite temperature, however, superfluidity in this model is 
destroyed by fermionic excitations, namely the breakup 



of pairs. The critical temperature is therefore of the same 
order as the pairing gap at zero temperature, consistent 
with the well known BCS relation 2Aq/T c = 3.52 in weak 
coupling. Clearly, such a picture is appropriate for weak 
coupling, where the transition to superfluidity is driven 
by the gain in potential energy associated with pair for- 
mation. By contrast, for sufficiently strong interactions, 
the supcrfluid to normal transition is instead driven by a 
gain in kinetic energy, associated with the condensation 
of pre-formed pairs rather than their thermal breakup. 
The critical temperature is then of the order of the de- 
generacy temperature of the gas and thus is completely 
unrelated to the pair binding energy. For a proper de- 
scription of the BCS-BEC crossover at finite tempera- 
ture and arbitrary coupling, we therefore need to go be- 
yond mean-field, including excitations, which drive the 
supcrfluid order parameter to zero without destroying 
the bound pairs altogether. On a formal level, this is ac- 
complished by the nontrivial wave-vector and frequency 
dependent term GT in the exact fermion self energy, as 
given in ([2730]) . 

For the numerical calculation we decompose the ther- 
modynamic potentials into a mean-field part and a cor- 
rection term according to fi = fix + Af2, S — Si + AS, 
etc.. The mean- field contributions have been derived in 
the previous subsection. While in these contributions the 
Matsubara-frequency sums have been performed explic- 
itly, the integrals over the wave vector remain and must 
be evaluated numerically. By subtracting the mean-field 
formulas ((2745)) - ((2747)) from the general formulas ((2~4"T)) - 
(|2.43p we obtain the correction for the grand thermody- 
namic potential 



An = -L d J ^^'{-^[^(^^^-^(^^^[^(^^^-^(k,^)-!]} 



(2.54) 



11 



the correction for the combination of the internal energy and the particle number 

At/ - ,iAN = - L d J l^Tr{[Gi(k,i J , > )- 1 +i/ii.,J[G(k,i.„)-G 1 (k,i l .„)]} 

and the correction for the entropy 

AS =/3L d J -0^ i5^TS:{-ln[Gf 1 (k,w n )- 1 G(k,£«; B )] + (-ifej n )[G(k,w n ) - Gi(k )Wn )]} 

-^/3i d / ^ ^Tr{-ln[r(K,fi„)/. 9o ] + [r(K,a i )/.g -l]} . 



(2.55) 



(2.56) 



In formulas (|2.54[) - (|2.56p the sums over the fermionic Matsubara frequencies uj n converge so that the regularization of 
Appendix [A] is not needed. However, the sums over the bosonic Matsubara frequencies Q, n are not well defined and 
must be regularized. Thus, for a numerical evaluation the formulas (|2.54p - (|2.56[) must be transformed further, which 
will be done in the next subsection. 



H. Renormalization of the thermodynamic potentials 

Since the interaction has zero range, the interaction strength g must be renormalized and replaced by the scat- 
tering amplitude g according to (|2.4| . In a first step we renormalize the mean-held formulas of the thermodynamic 
potentials. In (|2.48p - (|2.50[) the interaction strength go does not occur explicitly. The integrals are thus finite and no 
renormalization is needed for these formulas. However, the condensation energy (|2 . 5 1 [) contains two infinite terms, a 
divergent integral and the last term with the infinite factor 1/go, which compensate each other. By renormalizing the 
interaction strength we obtain 



-2L C 



d d k 1 



{2n) d 2 



IAI 



2e t 



A 



2 



(E k — (J,) — (e k -IJ.)-hr- ~ L — ( 2 - 57 ) 
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where both the integral and the last term are now separately finite. Note that the wave vector integrals in (|2.48|) - (|2.50j) 
and in (|2.57[) are finite in any spatial dimension d with 2 < d < 4. 

In a second step we renormalize the correction formulas. In the correction of the grand thermodynamic potential 
(T2.54[) we decompose ln[r(K, fl n )/go] = ln[T(K, Cl n )/g] + ln[g/ga\. The separated term ln[g/go] can be neglected 
because it does not depend on the Matsubara frequencies il n - Following the arguments of Appendix [5] the Matsubara- 
frequency sum of this term is zero. Thus, for the correction of the grand thermodynamic potential we obtain the 
renormalized formula 



An = -L d J ^ i^Tr{-ln[G 1 (k, W „)- 1 G(k,u;„)] + [G 1 (k, W „)- 1 G'(k,c l ; n )-l]} 
+ ^ d /|^^ Tr{ " ln[r(K ' ai)/5]} - 



Both terms of this formula are now finite. However, Eq. (|A4[) is needed to evaluate the second term. 

In correction (|2 . 55[) the second term must be renormalized. This can be achieved by the following sequence of 
equations 



Ld j0f-pY, n ^/go 1] } = L d J |* i £ Tr{r(K, O n ) x( K, fi n )} 

Ld J (2^ ^E Tr {^( k .^)G(k,c„)]} =L d J I^Tr{G 1 (k,o; n )- 1 [G(k ! a;„)-G 1 (k ) a;„)]} 



First, by using the Bcthc-Salpctcr equation (|2.33|) we write the integrand as a product of the matrix vertex function 
r and the matrix pair propagator x- Secondly, we express x m terms of the matrix Green function G by (|2.18[) . 
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interchange the orders of the integrals and sums, and combine T with one of the G into the self energy £ by (|2.40|) . 
Finally, we replace £ in favor of G and G\ by using the Dyson equation (|2.39p . The bosonic integral and sum are 
transformed into a fermionic integral and sum. Hence, the second term of (|2.55j) can be combined with the first term. 
As a result we finally obtain 

AU-fiAN=-^L d J i^Tr{[G 1 (k,a;„)- 1 +2^ n ][G(k, W „)-G 1 (k, W „)]} . (2.60) 

By considering (|2.38|) we explicitly prove 

Gi(k,w„) _1 +2ihu; n = G 1 {k,-(j n )- 1 . (2.61) 

Consequently, for the correction of the combination of the internal energy and the particle number we obtain the 
compact formula 

AU - fiAN = - l -L d J J^L i^Tr{G 1 (k,-c„)- 1 [G(k,c„)-G 1 (k, w , i )]} (2.62) 

which is essential for a stable numerical evaluation of the correction term. The Matsubara-frequency sum is evaluated 
by using (|A3[) . The wave- vector integral is finite. 

The correction of the entropy (|2.56|) is rcnormalized in an analogous way. Alternatively, we use the thermodynamic 
relation (|2.44|) . As a result we obtain 

AS =/3 L d J A^L 1 £ Tr {- ln[Gx(k, u^GQt, w„)] + [G x (k, cj n )- x G(k, u n ) - 1]} 

-\j Ld \ ^Tr{G 1 (k,-uj n )- 1 [G(k > w n )-G 1 (k,u; n )}} (2.63) 

-S^/^^^-HTCK.nnVff]}. 



The final results are the mean-field formulas (|2.48|) - 
(|2.50[) together with (|2 . 57[) and the correction formulas 
(|2~55)) . ([2l)2"]) . and pT55]) . In these formulas each term by 
itself is finite. Eventually, the thermodynamic potentials 
are obtained by adding the terms together according to 
ft = fix + Ail, S = Si + AS, etc.. 

I. Symmetry breaking and Thouless criterion 

The interacting fcrmion system is invariant under the 
symmetry transformation 

^(r)^e a VaW , ^+(r) e" iA ^(r) (2.64) 

which is related to a global change of phase of the 
fcrmion fields by A. The superfluid phase breaks this 
symmetry since the order parameter A is transformed 
as A — > e 2tX A. Clearly, however, the thermodynamic 
potentials must remain invariant under a global change 
of the phase both in the normal and in the superfluid 
state. In the superfluid, the free energy increase associ- 
ated with a slowly varying phase A(r) vanishes like (VA) 2 . 
By Goldstone' s theorem, this implies the existence of 
modes whose energy vanishes in the long wave-length 



I 

limit. For a neutral superfluid, this is the well known 
Bogoliubov- Anderson mode. It has a sound like disper- 
sion (j(k) = ck and is physically related to fluctuations 
of the phase of the order parameter. 

In technical terms, the existence of zero-energy collec- 
tive modes can be derived from Ward identities related to 
the symmetry transformation. By considering the grand 
thermodynamic potential n[G], in Ref.[3l|the Ward iden- 
tity 

Y, r xx'x>Y^ YY '=0 (2.65) 

YY' 

has been derived (see (2.57) in Ref. [3l|). Here 5\Y,xx> 
is the variation of the self energy under the transforma- 
tion (|2.64p with an infinitesimal phase change SX. This 
quantity may be interpreted as the generalized order pa- 
rameter of the system. On the other hand, T^-/ Y y> 1S 
the inverse vertex function. For a short-hand notation 
the indices X , X' and Y, Y' are used, which represent a 
combination X = (r, r, a, a) of the space variable r, the 
imaginary time r, the spin index a, and the Nambu index 
a. According to (|2.65j) the inverse vertex function 
may be interpreted as a linear operator which acts on the 
order parameter 5\T.. In the superfluid state the order 
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parameter is nonzero so that the inverse vertex function 
must have a zero eigenvalue, which is related to a zero 
energy collective mode. For superfluid fermion systems 
this fact is known as the Thouless criterion [39| . 

The Ward identity (|2.65|) has been derived for the 
exact theory. However, our present crossover theory is 
an approximation, based on a certain truncation of the 
exact functional which enters either into the Luttinger- 
Ward or the DeDominicis-Martin formalism. In general, 
such a truncated functional will not obey the Ward iden- 
tity. Indeed, we find that our inverse vertex function 
r~ , (K, O n ) obeys instead the equation 

£r-i,(K = O,n„ = 0)A Q , = C(|A| 3 ) (2.66) 

a' 

where (A a ) = (A, A*) (see (3.56) in Ref. HH). Taking 
the longitudinal part, this equation correctly describes 
the smooth evolution from a Ginzburg-Landau type de- 
scription of weak coupling BCS supcrfluids to a Gross- 
Pitaevskii like theory of a dilute, repulsive Bose gas [3lj . 
The transverse part, however, also gives a finite value 
on the right hand side of (|2 . 66[) in the limit K = and 
D, n = 0, thus violating the Ward identity by terms of 
order | A| 3 . As a result, the Thouless criterion is violated 
and there is no proper Bogoliubov- Anderson mode in our 
approach without a further modification (see below). 

Unfortunately, the violation of the Goldstonc theorem 
for continuous symmetries is a general property of con- 
serving approximations based on the Luttinger-Ward for- 
malism. This problem has been known for a long time 
for superfluid Bose systems [!o| and is sometimes referred 
to as the 'conserving-gapless dichotomy' [4l|, |42j in the 
literature. For the exact theory, a Ward identity holds 
for the inverse matrix boson Green function Gb, which 
reads 

^G s ; QQ ,(K = O,ft„ = 0)* B , tt <=0. (2.67) 

In the superfluid state the inverse matrix boson Green 
function has a vanishing eigenvalue. For superfluid bo- 
son systems this is known as the Hugenholtz-Pines the- 
orem [43l | . Conserving approximations, however, violate 
the Hugenholtz-Pines theorem. For example, this is true 
already for the lowest approximation, the well known 
Hartree-Fock-Bogoliubov theory. 

In our fermion system for strong attractive interac- 
tions v = I/Hfcl ^> 1, the fermions are bound into pairs. 
These pairs form a Bose system with an effective repul- 
sive interaction which, for a dilute system, is described 
by the exact scattering length add ~ 0.60 a of the four 
particle problem associated with dimer-dimer scattering. 
In the strong coupling limit, therefore, our crossover the- 
ory for interacting fermions must converge to an effective 
theory of repulsively interacting bosons, where both the- 
ories are based on the Luttinger-Ward formalism. From 
the analytical arguments in Refs.l29l.l3ll and also from our 
numerical calculations, we find that the crossover theory 



converges to the Hartree-Fock-Bogoliubov theory quickly 
for interactions v = 1 jkpa > 2. The boson order parame- 
ter and the matrix boson Green function Gs(K, f2„) 
can be identified with the order parameter A and the 
vertex function T(K, Q n ) according to [29|, [3l| 

* B = ±i[8Trela 3 ]- 1/2 A , (2.68) 
G B (K,n n ) = -[STrefa 3 ]- 1 ^^) . (2.69) 

The validity of these relations requires both strong cou- 
pling, but also low frequencies and momenta. Indeed, 
it is only at low energies where the composite particles 
behave like bosons. At higher frequencies or momenta, 
the composite nature of the pairs becomes visible. This 
becomes evident, for instance, in the different behaviour 
Gb ~ of a Bose Green function at large frequencies 
compared to that of the vertex function, which behaves 

— 1/2 

like r ~ Q n as a result of the two particle continuum 
associated with broken fermion pairs. Clearly, at large 
coupling constants » > 1, this continuum moves up to 
very large frequencies of the order of the binding energy 
e b ~ v 2 . 

We conclude that the violation of the Thouless crite- 
rion in our crossover theory is related to the violation 
of the Hugenholtz-Pines theorem in the Hartree-Fock- 
Bogoliubov theory for bosons to which our Luttinger- 
Ward formulation of the fermionic many-body problem 
converges at large coupling. In the following section, it 
will be shown that this problem may be solved by an 
appropriate modification of the coupling constant. In 
this manner, a self-consistent formulation of the many- 
body problem is possible which obeys Goldstone's the- 
orem and thus provides a correct description of both 
fermionic and collective, bosonic excitations along the 
BCS-BEC crossover. 



J. Modified coupling and gapless 
Bogoliubov- Anderson mode 

In the following, our aim is to modify the theory in a 
way which is consistent with the Thouless criterion, giv- 
ing rise to a gapless Goldstone mode in the whole regime 
of coupling strengths. If we require the Thouless criterion 

^^ = 0,^ = 0)^ = (2.70) 

ot' 

a further equation will be added to the self-consistent 
equations for the Green and vertex functions in Subsec. 
Ill Dl However, then another equation must be discarded 
or a further parameter must be introduced. We find 
that the Bethe-Salpeter equation (|2.35[) and the order- 
parameter equation (|2 . 3T[) can not be satisfied together, 
if (|2.70p is required. For this reason we modify the the- 
ory by introducing a modified scattering amplitude <? m od> 
which is determined by the modified order-parameter 
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equation 



A 



fju 



0) 



(2.71) 



We have solved the self-consistent equations together 
with the Thouless criterion (|2.70[) and the modified order- 
parameter equation (|2.71[) . The numerical effort is much 
less for the modified theory than for the original the- 
ory. Since the scattering amplitude g is related to the 
scattering length a, we obtain a modified dimensionless 
interaction strength u mo d = 1/fcFOmod- We find a dif- 
ference <5v mo d = v mo d ~ v in the range between 0.0 and 
-0.1. 

In order to obtain a consistent theory, we must check 
that the modification is compatible with the Luttinger- 
Ward formalism. We must find a modified grand ther- 
modynamic potential f2 mo d[G], so that the condition for 
stationarity (|2.9p yields the self-consistent equations with 
the modified order-parameter equation (|2. 71[) . For this 
purpose we consider the second term of (|2.16|) which 
reads 

n o [G] = L d 5o|^(O,0)| 2 = L d \A\ 2 /g . (2.72) 
We replace this term by the modified term 

tto,mod[G] = L d |A| 2 /5o,mod(|A|) (2.73) 

where 

A = 5o,mod(|A|)^(0,0) (2.74) 

is the modified order-parameter equation. The modified 
interaction strengths £/ ,mod = £/o,mod(|A|) and po.mod = 
5o,mod(|A|) depend on the order parameter |A|, are not 
equal, and differ from the bare interaction strength go- 
In order to apply the stationarity condition (|2.9[) we must 
consider the variation of (|2.73|) with respect to G. Since 
the modified parameter <7o.mod(|A|) depends implicitly on 
G via (|2.74[) , the chain rule of differential calculus must 
be applied. Eventually, the variation of (|2.73[) must have 
the form 



til 



'^[G]=L d J ^i^Trj^^k,^)} 
= L d [A ST(0, 0)* + A* 6f(0, 0)] . 



(2.75) 



By comparing the resulting terms with (|2.75|) , we obtain 
the differential equation 



d 



8\A\ 50,mod(|A|) 



2IAI 



d 



IAI 



9|A| 5 o,mod(|A|) 



(2.76) 



On the other hand, Eq. (|2.66|) implies that the Thouless 
criterion holds without modification if |A| = 0. Thus, we 
find 



5o,i 



5o for|A|=0 



(2.77) 



which is an initial condition for (|2 . 76|) . Eq. (|2.76|) can be 
integrated together with (|2.77|) . We obtain 



1 



1 



A| 2|A'|d|A'| 



50,raod(|A|) 50,mod(|A|) |A| 2 7 50.mod(|A' 

(2.78) 

The thermodynamic state of the interacting fermion sys- 
tem in the superfluid state is therefore determined by 
three parameters. We may choose the order parame- 
ter |A|, the chemical potential /i, and the interaction 
strength 170 for these parameters. Hence, the modi- 
fied interaction strengths go, mod = <?o,mod(|A|,/z, 50) and 
50, mod = 5o,mod(|A|,^,g ) are functions of these pa- 
rameters. While <7o,mod(| A|, /j,, go) is uniquely deter- 
mined by (|2.74p and the other self-consistent equations, 
<?o.mod(| A|, go) depends on the path in the parameter 
space when the integral (|2.78p is calculated. Since go and 
fi are external parameters of the theory, for a correct for- 
mulation of the modification these parameters must be 
kept constant. 

The modification is compatible also with the DcDo- 
minicis-Martin formalism. In this case the internal en- 
ergy U[G, T] includes the term (|2.72|) which must be mod- 
ified according to (|2.73p . The modification of the cou- 
pling constant go described by (|2.76|) - (|2.78|) is derived in 
an analogous way. 

Eqs. (|2~72]) - ((^751) describe the modification of the 
crossover theory in terms of the bare interaction parame- 
ters go, 5o,mod, and 50, mod- A renormalized version of the 
modification is obtained, if we replace the bare parame- 
ters by the renormalized scattering amplitudes g, g mo d, 
and g mod according to (pH]) . Eqs. (f2~76>(|2~T8l) are valid 
also for the renormalized scattering amplitudes without 
changes. From (|2 . T8[) we obtain 



1 



1 



|A| 2|A'|rf|A'| 

flmod(|A'|) 



ffmod(|A|) g mod (\A\) |AP J g mod (\A'\) ' (2J9) 

The renormalized modified order-parameter equation is 
defined by ()2.7ip . In order to obtain the modified for- 
mulas of the renormalized thermodynamic potentials in 
Subsec. Ill HI only a single change is needed. We must 
replace the energy constant (|2.57|) by 



O.mod 



= -2L 



d d k 



{2n) d 
L d |A| 2 (g 



2 (%-X 
mod 



(E k - n)- (e k 
2.9- 1 ) ■ 



(2.80) 



The other formulas ([2~58]) . ([2^62]) . and ([2~63]) remain un- 
changed. Since the renormalized scattering amplitude 
g is related to the dimensionless interaction parameter 
v = l/fci?a, we can transform (|2.79p into a dimensionless 
form. For Sv mod = v mod - v and <5"D mod = u mod - v we 
obtain 

(\A\/e F ,v) = 2<5?; m od(|A|/e F ,D) 

AA\/EF 

- (| A\/e F )- 2 / 6v mod {X, v) 2X dX . 
Jo 

(2.81) 
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FIG. 3: (Color online) Modifications of the dimensionless in- 
teraction parameter: red solid curve shows <5v m od and blue 
dashed curve depicts 5v mo d as a function of v for T = and 
|A| = |A |. 

While Sv mo( i is obtained directly from (|2.71[) by solving 
the self-consistent equations, <5w mo d is obtained by eval- 
uating the integral in (|2.8ip numerically. As a result we 
obtain modifications which are restricted to the interval 

- 0.1 < Sv mod < Sv mod < for |A| > . (2.82) 

In Fig. [3] the modifications <5u mo d and 6v mo d are shown 
as red solid curve and blue dashed curve, respectively, 
for T = and |A| = |A |. Clearly, the modifications are 
largest in the crossover region close to the unitarity point. 
At finite temperature for increasing T the order parame- 
ter |A| and the modifications 5v mo d and <5vmod decrease 
together. Eventually, for |A| = the modifications are 

femod = <5'5mod = 0. 

In the previous subsection we have argued that for 
strong attractive interactions the fermions are bound into 
pairs. Our crossover theory for the interacting fermion 
system then converges to a Luttinger-Ward type the- 
ory for interacting bosons. Since the modified version 
of our theory obeys the Ward identity, its strong cou- 
pling limit necessarily leads to a description of dilute, 
repulsive bosons which has the correct linear spectrum 
of excitations at low energies. It turns out, that the lim- 
iting theory here is the Luttinger-Ward version of the so- 
called Shohno theory [3, HH which is equivalent to the 
more well known Popov approximation. While we have 
not been able to derive the Shohno-Popov theory analyt- 
ically from the Luttinger-Ward functional of the original 
fermionic model, we find a quick convergence numeri- 
cally in all thermodynamic quantities for dimensionless 
couplings v > 2. Considering the entropy, in particu- 
lar, the Shohno-Popov theory gives rise to the standard 
expression 

S = L d J |^{(l +n £)ln(l+n£)-n|lnn£} (2.83) 



for the entropy of a nonintcracting gas of bosonic quasi- 
particlcs with the standard distribution function 

n£ = l/[eMP[Ei - Ms]) - 1] ■ (2-84) 
The corresponding spectrum of excitation energies 

Ei-fiB = [(h 2 K 2 /2m B ) 2 + (h 2 K 2 /2m B )2g B \y B \ 2 ] 1/2 

(2.85) 

has the well known form of a Bogoliubov spectrum with 
a temperature dependent condensate density n B fi = 
\^ B \ 2 and a positive Bose-Bose scattering amplitude g B . 
Within our approximation, we have g B = 2g. i.e. the ex- 
act dimcr-dimcr scattering length add ~ 0.60 a is replaced 
by its Born approximation result a d ^ = 2a [46[. The ef- 
fective mass and chemical potential take their obvious 
values m B = 2m and \i B = 2fi + eb where = H 2 /ma 2 is 
the two-particle binding energy on the BEC-side of the 
crossover, where a > 0. The order parameter is given 
by (|2.68[) . Other thermodynamic quantities are obtained 
using more complicated expressions, which are not pre- 
sented here in detail [47j . As will be shown explicitely 
in the following sections, the numerical results for the 
critical temperature or the entropy converge quickly to 
that of the Shohno-Popov theory for coupling strengths 
v > 2. 



III. NUMERICAL RESULTS 

Following the detailed discussion of the formalism used 
to describe the thermodynamics of attractively interact- 
ing fermions at arbitrary coupling and temperature, we 
now present numerical results which cover both the nor- 
mal and superfluid regime. These results require a so- 
lution of the self-consistent equations determining the 
Green function G and the vertex function T, which are 
scalars above, and two-by-two matrices below the crit- 
ical temperature. An iteration procedure is performed 
where a numerical Fourier transformation is needed to 
transform the functions back and forth. Since the Green 
function G, the vertex function T, and the related func- 
tions E and M are singular at small values of r and r 
and also exhibit significant variation over several orders 
of magnitude, the numerical Fourier transformation is 
quite challenging. In practice, the variables need to be 
discretized on logarithmic scales. Standard procedures 
like the fast Fourier transformation are therefore not ap- 
plicable. The basic principles of our special numerical 
Fourier transformation are described in Appendix [B] 

A. Critical temperature 

The crucial quantity which determines the overall 
structure of the phase diagram is of course the critical 
temperature T c for the transition to a superfluid. This 
temperature is known analytically only in the extreme 
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BCS- and BEC-limit. In the BCS-limit k F \a\ -> 0, where 
the average distance between the fermions is much larger 
than the magnitude of the scattering length, the standard 
solution of the gap-equation for an attractive pseudo- 
potential gives a critical temperature 



=7e 



t (bcs) = - £F exp(-7r/2fc F |a|) 



(3.1) 



with 7£ = 0.5772 . . . Euler's constant. r £^ BCS ^> j s expo- 
nentially small on the characteristic scale of the Fermi 
energy. Since typical Fermi temperatures in cold gases 
are of the order of micro-Kelvin, the BCS-rcgime is in 
practice hardly attainable in these systems. 

The leading order corrections to the BCS-result in 
an expansion in the small parameter kp\a\ <C 1 have 
been determined a long time ago by Gorkov and Melik- 
Barkhudarov (48|. They arise from induced interactions, 
where one fermion sees the polarization in the Fermi 
gas due to a second fermion. The density induced in- 
teraction changes the dimensionless coupling constant 
N(0)g = 2k F a/n of the BCS-theory to [|| to 



g^g + g 2 N(0) 



1 + 2 In 2 



(3.2) 



where N(0) = mkp/2Tr 2 h 2 is the standard density of 
states per spin at the Fermi energy. Since the addi- 
tional contribution to the two-body scattering amplitude 
g < is positive, the induced interactions weaken the 
attractive interaction between two fermions in vacuum 
and lead to a reduction of the transition temperature by 
a factor (4e) _1 / 3 ss 0.45. The nonanalytic dependence 
of the BCS-transition temperature on the dimensionless 
coupling constant kpa thus give rise to a finite change in 
the prefactor in ((33]) from the BCS value 0.61 to 0.28, 
even though the contribution of induced interactions is 
of order kpa compared to the bare interaction. 

On the BEC-side, the zeroth order result for the critical 
temperature is obtained from the value 



2/3 



t (bec) = 3 31 §_ = o.218 e F 



TUB 



(3.3) 



obtained for an ideal Bose gas with density ns — n/2 and 
mass rns = 2m. The leading corrections to this result 
arise from the residual interactions between the stro ngly 
bound bosonic dimers. As shown by Petrov et al. [2fl[50|. 
these interactions can be described by a positive dimer- 
dimer scattering length add ~ 0.60 a. With the quite 
plausible assumption, that the total potential energy in 
a dilute gas of dimers is the sum of its two-body interac- 
tions, the scattering length of the four fermion problem 
determines the corresponding interaction constant in the 
theory of a weakly interacting Bose gas in the regime of 

1 /3 

a small gas parameter n B add <C 1, where Bogoliubov 
theory is applicable. The exact dependence of the criti- 
cal temperature of the dilute, repulsive Bose gas on the 
interaction strength has been calculated only in recent 



years. To lowest order in the interaction, the shift is 
positive and linear in the scattering length [5lj . 



Tc/T c 



(BEC) _ 



1 



1/3 

cn R a d d 



(3.4) 



with a numerical constant c ~ 1.31 [52|, [53[. As a result, 
the evolution of the critical temperature in the homoge- 
neous case as a function of the dimensionless coupling 
constant v = 1 /kpa necessarily exhibits a maximum, 
since the asymptotic ideal Bose gas result is approached 
from above. Such a maximum has been found in the 
early calculations of T c along the BCS-BEC crossover by 
Nozieres and Schmitt-Rink [54| and by Randeria et al. 
[55| . The precise height and location of this maximum, 
however, has not been determined so far in a quantita- 
tively reliable manner. Given that our present theory 
exhibits a first order transition, there is a range of multi- 
valuedness of the thermodynamic potentials as a function 
of temperature. This regime is bounded in Fig. [5] by the 
upper and lower T c curves respectively. The lower T c 
curve (shown as the red dashed line) which is monotonic 
in v coincides with the T c curve previously calculated 
[30j | by implementing the Thouless criterion coming from 
the normal fluid side. In a situation, where a true first 
order transition is expected, we would need to perform 
a Maxwell construction to obtain the proper transition 
line. As was discussed above, however, the first order 
transition is an artefact of the approximations involved. 
In particular the spectrum of excitations right at T c is free 
particle like in our approximation rather than ujk ~ K 3 ^ 2 

In order to determine the proper critical temperature 
within our approximation , we have used two essentially 
equivalent criteria: the fact that the exact entropy is con- 
tinuous at T c suggests that our best approximation for 
the critical temperature is where the jump in the entropy 
between the two branches characterising the superfiuid 
and the normal regime has a minimum. Essentially the 
same value is obtained by defining T c through the cri- 
terion that it is the maximum temperature at which the 
order parameter A(T) is nonzero. Remarkably, these cri- 
teria lead to a critical temperature (shown as 0i upper) in 
Fig. 2]) which exhibits a maximum on the BEC-side of the 
crossover around v « 1 as expected on general grounds. 
Moreover, our theory predicts the correct asymptotic 
functional form (|3.4[) of the T c -enhancement in the BEC 
limit v 3> 1. Even though the dimer-dimer scattering 

( B) 

length a dd = 2 a and the prefactor c sa 0.58 of our ap- 
proximate Popov-type theory differ from the exact values 
Odd = 0.60 a and c w 1.31, respectively, the agreement of 
our theory with the exact result is very good (see Fig. QJ. 

A quite sensitive test of the quantitative reliability of 
our present result for the critical temperature at arbi- 
trary coupling is provided by a comparison with the re- 
cent, rather precise numerical results right at the unitar- 
ity point by Burovski et al. [28[. In fact, our result for the 
dimensionless ratio T c /sp ps 0.16 of the critical temper- 
ature in units of the bare Fermi energy, which is one of 
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FIG. 4: (Color online) 6i lower) (red dashed line) and 8 { c upper) 
(solid black line, identified as T c ) compared with the Shohno 
result (blue dotted-dashed line) with = 2 a and the 

exact (QMC) result (light-blue squares) with AT c /Tbec — 

1/3 

«!j add and c = 1.31 and add = 0.60 a. Yellow dashed line 
and green triangles show the BCS result without and with 
Gorkov and Melik-Barkhudarov corrections. 
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FIG. 5: (Color online) S(T) at various interaction strengths 
v. 



the universal numbers of the BCS-BEC crossover prob- 
lem (see Subsec. UlI Bl below), agrees precisely with the 
numerical results of Burovski et al. within the given er- 
ror bars. As will be shown below, a similar rather precise 
agreement is obtained with other thermodynamic quanti- 
ties, except for the chemical potential. Thus, even in the 
absence of a small parameter which would allow to con- 
trol our theory systematically in the crossover regime, the 
agreement with the numerical results at unitarity gives 
us confidence that the approach is quantitatively reliable 
at arbitrary coupling strengths. 

In Fig. [5] the temperature evolution of the entropy is 
shown for various coupling parameters v. Here the mul- 
tivalued character is clearly seen which reflects the first- 
order transition. Furthermore, three-dimensional plots 
of the entropy and of the pressure are presented in Figs. 
|5]and respectively. In both figures a rather sharp drop 




FIG. 6: Entropy as a function of 6 and v obtained using (|2.50[) 
and (f2763|) . 
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FIG. 7: Pressure as a function of 6 and v obtained using (|2.48p 
and (p38l) . 



is observed in the crossover region from weak coupling 
v < —1 (fcrmionic regime) to strong coupling v > +1 
(bosonic regime) . In the weak coupling limit v <C — 1 the 
results of the nearly ideal Fermi gas are approached which 
are defined by the BCS formulas <|2Tl8]) - (f2T50|) and f237|) . 
On the other hand, in the strong coupling limit v +1 
the results of Shohno's mean-field theory are approached. 
While the strong-coupling entropy is defined by (12.831) . 
the other thermodynamic quantities are defined by more 
complicated formulae |47j |. For v > 1.0 and very low 
temperatures the pressure is nearly zero which reflects a 
special property of Shohno's mean-field theory of weakly 
interacting bosons. At high temperatures T 3> £_f the 
entropy, the pressure, and the related thermodynamic 
quantities approach the Boltzmann limit. 

Fig. [5] shows the order parameter which vanishes ex- 
ponentially A(T = 0)/sf — ► (8/e 2 ) cxp(7ru/2) according 




FIG. 8: (Color online) 3d view of the order parameter. 



to the well known BCS result for v <C —1.0. In the op- 
posite limit of strong coupling the behaviour can be de- 
rived from [in — ► —A 2 /2g which reflects the fact that the 
fermion chemical potential in the strong coupling limit 
is governed by the potential (i.e. binding) energy. This 
yields A(T = Q)/ep — > y/ (16/37r)w with the square root 
behaviour clearly visible in Fig. [5] Near T c the gap func- 
tion displays the multivalued behaviour characteristic of 
a first-order transition. 

At low temperatures the entropy has to vanish, in ac- 
cordance with the third law of thermodynamics. The 
way it does, is in fact universal along the full BCS- 
BEC crossover. Indeed, at low temperatures, the two- 
component Fermi gas is in a superfluid state, independent 
of the strength of the attractive interaction. On quite 
general grounds therefore, the low lying excitations above 
the ground state are sound modes of the Bogoliubov- 
Andcrson type. These modes give rise to an entropy 



(3.5) 



which vanishes like T 3 for arbitrary coupling strength. 
The associated sound velocity c is constant at low T and 
may be determined from the pressure via mc 2 = dp/dn. 
Fig. [5] displays (c/vp) 2 at T = as a function of cou- 
pling strength with vp the Fermi velocity. The dilute 
interacting Fermi gas limit (c/vf) 2 = (l + 2/(nv))/3 and 
the BEC limit (c/vp) 2 = kpcidd/ '(6V) for add = 0.60 a are 
represented by the blue squares and the green triangles 
respectively. The red triangles are obtained by extend- 
ing the expression of the ground state energy of a dilute 
weakly interacting Fermi gas [5?], HI] with the help of a 
Pade approximation to the strong coupling regime [1, [g[ 
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FIG. 9: (Color online) Isothermal sound speed mc 2 = dp/dn 
as a function of u = 1/fcpa for T = 0. The different curves 
are explained in the main text. 



and the thermodynamic identity 
m dn 



n2 dEjN^ 



dn 



(3.7) 



Obviously the present crossover theory provides a very 
good description of the equation of state and sound veloc- 
ity except in the regime v > 1, where our results under- 
estimate both the pressure and its density dependence. 

In principle we should be able to independently obtain 
c from the low entropy asymptotics (|3.5[) . Our numerical 
results are consistent with S(T) ~ T 3 , however they are 
not precise enough at such low temperatures, to extract 
the sound velocity in this manner. 



B. Thermodynamics in the unitarity limit 

After presenting the results for the critical temperature 
and the thermodynamics at arbitrary coupling, we now 
turn to a more detailed discussion of the unitarity limit, 
where the scattering length is infinite. This particular 
line in the phase diagram has received a lot of atten- 
tion recently. In particular, precise numerical results are 
available at this point [28| , which provide a sensitive test 
of analytical approaches to the crossover problem. 

As has been mentioned before, the Fermi gas at infi- 
nite scattering length v — is rather special since the 
only relevant length and energy scales remaining in the 
problem are the Fermi wave length set by the density and 
the Fermi energy e p , provided we remain within the zero 
range pseudopotential approximation. The free energy 
has a simple scaling form 



F(T,V,N) = f(8)Ne F 



(3.8) 



In particular, there are a number of universal ratios which 
characterize the crossover problem right at the unitarity 
point, both at zero temperature and at T c . Examples, 
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FIG. 10: Internal energy at unitarity as a function of temper- 
ature calculated using p,49[) and (|2.62[) . The dashed curve is 
obtained from the calculated pressure using the scaling for- 
mula U = \pV valid at unitarity. 
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FIG. 11: The single particle chemical potential at unitarity 
as a function of temperature obtained from the number con- 
servation constraint (|2.26p . 



which will be determined below, are the chemical poten- 
tial and the internal energy in units of the Fermi energy 
or the entropy per particle at T c . In addition, also the 
gap for single particle excitations or the condensate frac- 
tion at zero temperature are universal at the unitarity 
point. 

Fig. llOl shows the temperature dependence of the inter- 
nal energy calculated in two different ways. The solid line 
is our numerical result which is compared with the inter- 
nal energy (depicted as the dashed line) as obtained from 
the numerically calculated pressure p = — Q/V via the 
scaling relation U = 3pV/2 valid at the unitarity point. 
Our numerical results display perfect scaling above T c . 
The scaling violation below T c is a consequence of the 
modification of the theory. In order to preserve the con- 
serving nature of our theory while obeying the Thouless 
criterion an extra length scale a mo d had to be introduced 
leading to a modified dimensionless interaction strength 
fmod = 1/fcFamod with Sv mo(i = v mod -v in the range be- 
tween 0.0 and —0.1 with u mo d 7^ for v = (see Subsecs. 
|nTland[irj]for details). 

Fig. [TT] displays the behaviour of the chemical potential 
n(T) as a function of temperature. Using fi(T) in a local 
density approximation 



ti = ti h [n(r),T(r)/T F ] +V(r) 



(3.9) 



with fXh the chemical potential of the homogenous case we 
can calculate the density profiles of harmonically trapped 
ultracold gases at unitarity [59| . We have also checked 
the convergence of our fJ-(T) to the hi gh t emperature ex- 
pansion obtained by Ho and Mueller [60| which however 
only occurs for T 3> ef> 

Note that below T c the chemical potential n(T) is an 
increasing function of T. This perhaps counterintuitive 
result can be understood quite easily from the fact, that 
the low temperature thermodynamics is determined by 




FIG. 12: Entropy at unitarity as a function of temperature. 



the Bogoliubov- Anderson mode. As argued in the previ- 
ous section, this leads to an entropy which vanishes with 
a power law ~ T 3 . Fig. [12] displays the entropy at unitar- 
ity as a function of temperature. Now, at a given volume, 
there is a Maxwell relation of the form 



dfi 

df 



N,V 



dS 
6W 



T,V 



(3.10) 



which connects the temperature dependence of the chem- 
ical potential to the density dependence of the entropy. 
Using the universal result (|3.5[) for the low temperature 
entropy, this relation shows that at low temperatures the 
chemical potential exhibits a T 4 dependence with a pref- 
actor determined by 



df 



3S d 2 p 



> 



(3.11) 



2Vmc 2 s dn 2 

Obviously, this argument is not confined to the unitarity 
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O.lr 



0.08 -. 



0.06 



H^O.04- 



0.02 



T IT 

BEC v=2 



TABLE II: Comparison with diagrammatic determinant Mon- 
te Carlo (Burovski et al. (2£|), quantum Monte Carlo (Bulgac 
et al. |27fl), e = 4 — d expansion (Nishida and Son [2J, |2(|), 
Borel-Pade approximation connecting an expansion in e = 
4 — d and one in s — d — 2 [2(|) and a 1/iV expansion (Nikolic 
and Sachdev)Gl| at T = T c . 







T c /e F 




U/Ne F 


P/uef 


S/N 




Bulgac 


0.23(2) 


0.45 


0.41 


0.27 


0.99 




Burovski 


0.152(7) 


0.493(14) 


0.31(1) 


0.207(7) 


0.16(2) 




Nikolic (N = l) 


0.136 


0.585 


0.164 


0.109 






Nishida ( e = l) 


0.249 


0.18 


0.212 


0.135 


0.698 




Borel-Pade 


0.183 


0.294 


0.270 


0.172 


0.642 


4 


present work 


0.160 


0.394 


0.304 


0.204 


0.71 



FIG. 13: Temperature reduction on performing an isentropic 
sweep across v = from v = 2 to v = —2. 



TABLE III: Comparison with fixed node Green function 
Monte Carlo (Astrakharchik et al. [13] and Carlson et al. [l6( |) 
at T = 



TABLE I: Recent experimental results for (3 compared with 
calculated values. 











Bartenstein et al. [61] 


-o.esi^ 


Experimental 


Bourdel (2004) et al. [62] 


-0.64(15) 


results 


Duke (2005) [63J 


-0.49(4) 




Partridge et al. [64j 


-0.54(5) 




Astrakharchik et al. [17] 


-0.58(1) 




Carlson et al. [161 


-0.56(1) 


Calculated 


Hu/Liu/Drummond [67] 


-0.599 


values 


Perali et al. [65] 


-0.545 




Pade approximation [8, 9} 


-0.67 




present work 


-0.64 



point, showing that the chemical potential at low T has 
a behaviour fJ,(T) = /Lt(0) + 0{T A ) for arbitrary coupling 
strengths along the BCS-BEC crossover. A well docu- 
mented quantity which determines the density profile of 
dilute fcrmions in a trap at unitarity and T = is the so 
called P parameter defined via 



n(T = O)=e F (l + 0) 



(3.12) 



Our value of (3 ~ —0.640 is very close to j3 = —0.67 ob- 
tained via simply Pade approximating the weak coupling 
result for the ground state energy |§j9j and the experi- 
mental results of Bartenstein et al. 6l[ /3 = —0.68^'^ 
and Bourdel et al. [62| (3 = —0.64 ±0.15 but smaller than 
the results obtained at Duke [63[ , at Rice [64| and recent 
QMC results [H, Gjj (see Table HJ). Evidently, there is 
still considerable uncertainty in both the experimental 
and theoretical results. 

A promising route in the direction of thermometry for 
trapped gases is provided via the reversible adiabatic 
(isentropic) sweeps [1, [6l[ from the BEC limit. In Fig. 
[T3]we depict the resulting changes in temperature when 
moving across the unitarity limit for the homogenous 







U/Ne F 


P/ne F 


A/e F 


Astrakharchik 


0.41(2) 


0.25(1) 


0.17(1) 




Carlson 


0.43(1) 


0.26(1) 


0.17(1) 


0.54 


present work 


0.36 


0.21 


0.15 


0.46 



case. For the trapped case this cooling mechanism was 
first advocated by Carr et al. [32j and recently quantita- 
tively refined by Hu et al. [66[ ■ Finally to facilitate quan- 
titative comparison with various Quantum Monte Carlo 
results we have collected available data from the litera- 
ture presented in Table U and |TTJ Apart from the value 
for T c which is explicitely quoted in the paper by Bulgac 
(with errors) we have estimated the remaining quantities 
from their presented results and utilized scaling to fill in 
the missing data below. 

The T = results are fixed node QM results by As- 
trakharchik et al. and Carlson et al. [l(|. Note that 
our result for A/ep is close to the value Aqmb/£_f = 
(2/e) 7 / 3 = 0.49 obtained by a naive extrapolation of the 
Gorkov Mclik-Barkudarov result to fc^a = 00. 

At T c our results are in very good agreement with those 
of Burovski et al. except for the value of the dimension- 
less chemical potential (jl/sf and that of the entropy per 
particle at T c . Now Burovski et al. have obtained their 
values for the pressure p/tief and the entropy S/N indi- 
rectly from the internal energy and the chemical potential 
by using 3pV = 2U right at unitarity and the Gibbs- 
Duhem relation. The different results for the chemical 
potential then entail the considerable discrepancy in the 
value of S/N at T c . Within our numerical scheme, the 
chemical potential is the most directly - via (|2.26p - ob- 
tainable quantity among the thermodynamic data. In 
light of the excellent agreement of all other quantities 
with the numerical results of Burovski et al., the discrep- 
ancy for the chemical potential is thus quite surprising. 
Indeed, we believe that our values for both the chemi- 
cal potential and the entropy, for which the validity of 



the Gibbs-Duhcm relation and of 3pV = 2U at unitarity 
have been checked independently, are rather close to the 
exact results. This point of view is supported by consid- 
ering the evolution of the entropy per particle right at 
the critical temperature as a function of the dimension- 
less coupling. In the BCS limit, the entropy associated 
with single particle excitations can be calculated from the 
exactly soluble reduced BCS-Hamiltonian and is given by 
the standard mean-field expression (|2.50[) . At the critical 
temperature this entropy coincides with that of an ideal 
Fermi gas 

S(T C )/N= (ir 2 /2)(T c /T F ) . (3.13) 

Since the ratio T c /Tp is exponentially small in the 
weak coupling limit, the entropy (|3.13[) associated with 
fermionic excitations is dominant compared to the contri- 
bution arising from the collective Bogoliubov-Anderson 
mode. Indeed, extrapolating the corresponding low tem- 
perature entropy (|3.5[) associated with collective excita- 
tions up to the critical temperature gives rise to a contri- 
bution of order (T c /Tp) 3 , which is negligible compared 
to 1(313]) . 

At very large coupling strengths, the strongly bound 
fermion pairs form an eventually ideal Bosc gas, for which 
the entropy per particle right at T c can again be de- 
termined analytically. Recalling, that the number of 
bosons Nb = N/2 in this limit is just half the number of 
fermions, we obtain a universal number 

S(T C )/JV= ^|=0.6417.... (3.14) 

As is evident from Fig. [Ml where the complete evolution 
of the ratio S(T C )/N is shown as a function of the dimcn- 
sionlcss coupling parameter v, the limiting value of the 
ideal Bose gas is in fact not far from the entropy which is 
obtained from the Shohno-Popov theory of noninteract- 
ing bosonic quasiparticles in the range v > 1, according 
to ((2783]) . 

It is interesting to note, that the entropy per particle 
right at T c exhibits a maximum as a function of the cou- 
pling constant of order S(T C )/N f» 0.78 around the same 
coupling, where the critical temperature exhibits a max- 
imum. Considering the smooth evolution of S(T C )/N as 
a function of v, the value S(T C )/N w 0.16 at unitarity, 
which is deduced from the results of Burovski et ai, ap- 
pears to be far too small. On the other hand, the result 
S(T C )/N w 0.99 obtained by Bulgac et al. seems to be 
too high. 

IV. DISCUSSION AND CONCLUSION 

In conclusion let us summarize what has been achieved, 
mention shortcomings of the present approach and indi- 
cate possible future extensions. 

The formal basis of our results is a self-consistent, con- 
serving theory, which is based on an approach due to 




V = l/k_a 

F 



FIG. 14: Entropy at T c as a function of v = 1/kpa. Numer- 
ical result (solid) line obtained with (|2.5Up and (|2.63p com- 
pared with the limiting results: the BCS mean-field result 
(triangles) from (|2.50|l and (dashed line) from (|3.13|l and the 
Shohno-Popov result (dotted-dashed line) from (|2.84[) . 



Luttinger-Ward and DcDominicis-Martin, in which the 
exact one- or two-particle Green functions serve as an in- 
finite set of variational parameters. In order for this ap- 
proach to provide consistent thermodynamic results it is 
essential that the Green functions satisfy self-consistency 
conditions which reflect the stationarity of the appro- 
priate thermodynamic potentials. Approximate formula- 
tions, in which free Green functions are replaced by full 
ones according to a choice of Go Go, GGo or GG will in 
general not obey conservation laws or exact thermody- 
namic identities, in contrast to the $ derivable formu- 
lation presented here. The stationarity conditions were 
also crucial for the proof of thermodynamic equivalence 
of the Luttinger-Ward with the DeDominicis-Martin for- 
malism on the level of our approximate functional for 
the grand canonical potential or the entropy. In fact, to 
our knowledge, the theory presented here is the first con- 
crete application of the DcDominicis-Martin formulation 
to the fermionic many-body problem. 

An important point, we want to emphasize, is the nec- 
essarily self-consistent nature of the formalism. Indeed, 
within the Luttinger-Ward or the DeDominicis-Martin 
formulation an approximate functional for the grand 
canonical potential f2[G] or the entropy S[G, T] is made 
stationary by determining the space- and time-dependent 
Green and vertex functions from the variational condi- 
tions (|2.9p and (|2.13p respectively. The solution of these 
equations necessarily leads to a self-consistent mutual de- 
pendence of the various Green functions. Self consistency 
is thus reached precisely at the stationary point of these 
functionals. At this point, equations (|2.9p and (|2.13p are 
valid, conditions which arc necessary for the theory to 
give consistent thermodynamics, as pointed out e.g. in 
the context of equation (|2.27[) . 

A well known shortcoming of conserving approxima- 
tions is the dichotomy with the gaplcss nature of the 
collective modes, which reflects the broken continuous 
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symmetry of the superfluid state. For the present theory 
the formal reason for this dichotomy is a violation of the 
Ward identity resulting from the global gauge symmetry 
of the exact theory. In order to overcome this problem 
an extension of the theory was introduced which forces 
the gapless nature in the symmetry broken phase while 
remaining ^-derivable at the same time to maintain the 
conserving property. 

We have provided quantitative results for essentially 
all thermodynamic properties at temperatures below half 
the Fermi temperature, thus covering the relevant regime 
of the degenerate gas. Overall our results agree remark- 
ably well with recent numerical calculations at the uni- 
tarity point giving confidence that our approach is quan- 
titatively reliable over the full range of couplings between 
the BCS- and the BEC-limit. In particular, we provide 
concrete predictions for a number of universal ratios char- 
acterizing the unitary Fermi gas both at T = and at 
T = T C . 

The extensive numerical work entering the solution of 
the stationarity constraints and thermodynamic poten- 
tials is reflected most clearly in the three dimensional 
plots of the entropy Fig. pressure Fig. [7] and the or- 
der parameter Fig. [5J Most noteworthy are the quite 
abrupt change from fermionic to bosonic character for v 
in the interval — 1 < v < +1 which are most obvious in 
the entropy and pressure and the quick convergence to 
a Shohno-Popov theory of noninteracting bosonic quasi- 
particlcs for v > 1. 

An initially unexpected result, which is clearly visi- 
ble in the numerical data, is the fact that our superfluid 
phase transition is weakly first order, instead of being 
continuous as it should be. The origin of this failure to 
capture the critical behaviour correctly is found in the 
Shohno-Popov theory, which is obtained from our ap- 
proach in the limit d>1. The Shohno-Popov theory of 
a dilute, repulsive Bose gas generalizes the Bogoliubov 
theory to finite temperatures. It takes into account the 
thermal depletion of the condensate by including the ef- 
fect of bosons with finite momentum K in the particle 
number equation. Long ago Reatto and Straley [44| an- 
alyzed Shohno's theory in a self-consistent formulation 
and obtained a first-order superfluid transition. Phys- 
ically, the origin of the associated entropy jump is the 
collapse of the single particle spectrum right at the transi- 
tion. Indeed, within the Shohno-Popov theory, the single 
particle spectrum changes from initially linear to initially 
quadratic on raising the temperature through T c . As a 
result, the density of states is changed from a e 2 depen- 
dence below T c to the free particle y/e result right at and 
above T c . The associated drastic increase in the available 
phasespace leads to a jump in the entropy. 

For a purely bosonic system, a proper treatment of 
the behaviour near the critical point was recently given 
by Baym and coworkers [5l|. Baym and Holzmann [68j 
showed that a change of the spectrum for long wavelength 
excitations occurs right at T c . This hardening of the spec- 
trum (the low K spectrum is of the form K Q with a < 2) 



leads to the required reduction in the density of states to 
render the superfluid transition continuous. The subtle 
low K correlations necessary for this change in spectrum 
are clearly missing in our self-consistent approach. 

The BCS-BEC crossover being continuous however 
implies that the first order result also pertains to the 
V <C — 1 limit of our theory. We have checked that at the 
transition the discontinuities of all thermodynamic quan- 
tities are ~ exp(— C\v\) for v <C — 1 [13 • The associated 
difficulties of a proper treatment of bosonic excitations do 
not occur in the reduced BCS hamiltonian which neglects 
collective modes altogether, resulting in a continuous su- 
perfluid transition. To correctly account for the critical 
regime AT/T C — > our theory would need to be extended 
to treat the feedback between different bosonic modes ac- 
curately. Bickers and Scalapino have shown that this 
requires the incorporation of single particle self consis- 
tency and two particle self consistency on the same level 
of approximation. This may be achieved via so called 
parquet resummations. Currently, however, the inclusion 
of these contributions appears extremely challenging. A 
systematic and analytically accessible description of the 
crossover which is uniformly valid in both the normal and 
superfluid regime and which gives a proper account of the 
critical behaviour is provided by a 1 /./V-expansion as re- 
cently shown by Nikolic and Sachdev [l5[. This method 
can in fact be extended in a straigthforward manner to 
the case of unbalanced spin populations, a subject which 
has attracted a lot of attention very recently [U [7(| • 
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APPENDIX A: REGULARIZATION OF 
DIVERGENT MATSUBARA-FREQUENCY SUMS 

In our formulas of the thermodynamic potentials most 
sums over Matsubara frequencies arc not well defined. 
The functions which are summed do not decay to zero 
fast enough so that the Matsubara-frequency sums di- 
verge. However, this problem can be fixed. To do this 
we first perform a Fourier back transformation to obtain 
a function in terms of the imaginary time r. Then we 
take the limit r -> -0 or r -> +0 which is finite and well 
defined. 

We must distinguish between fermion and boson 
functions which have different Matsubara frequencies. 
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Fermion functions are of the type 



A(k,ui n ) 



A(k,u n ) 



B(k,w„) 

A{k,uj n y 



(Al) 



where A(k, u> n ) may be either ^4(k,w„) = — ln[G(k, w„)] 
or A(k,u n ) = [G (k,w„)- 1 G(k,w„) - 1]. (Note that 
the lower row of the matrix (|Alj) has the opposite sign 
than the lower row of the matrix (|2.17[) . The reason ist 
that in the terms of the thermodynamic potentials always 
an even number of fermion Green functions is multiplied 
together.) In this case we define 



\ ]T Tr{,4(k, ,„)} = 2 A(k, r = -0) 



(A2) 



where we assume that A(k, r) is real. Similarly consider 
a bosonic function of the form 



A(K O ^_ M(K,O n ) B(K,O n ) 



(A3) 



with A(K,0„) = r(K,Q n ) 

or A(K, il n ) = — ln[r(K, f2„)]. In this case we define 

i^Tr{A(K,cj n )} = 2.A(K,T = -0) , (A4) 







On 



where we assume that «4(K, r) is real. 

In some terms of our formulas the fermion function 
A(k, oj n ) or the boson function A(K, Sl n ) is proportional 
to the unit matrix 1. In this case the Fourier backtrans- 
form A(K, t) is 8f{t/K) or 8b{t/K), respectively. Hence, 
the related Matsubara-frequency sums (|A2|) or (|A4[) are 
zero. 



APPENDIX B: NUMERICAL FOURIER 
TRANSFORMATION 

The special numerical Fourier transformation has been 
invented long time ago by one of the authors in a differ- 
ent context in order to solve the mode-coupling equa- 
tion for the liquid-glass transition [7l| . In this case 
relaxation phenomena are considered on a logarithmic 
time scale over many decades, starting at microscopi- 
cally short times and extending up to very long macro- 
scopic times. Thus, a Fourier transformation is needed 
which can handle functions with features on logarithmic 
time and frequency scales extending over ten and more 
decades. Clearly, a standard fast Fourier transformation 
can not be applied because a constant step width would 
be needed. Rather the function to be transformed has 
been discrctized on a logarithmic scale and interpolated 
by cubic spline polynomials. Since for polynomial func- 
tions the Fourier integrals can be evaluated exactly, we 
end up with a transformation formula which depends on 
the spline coefficients of the function. 

Later this special numerical Fourier transformation has 
been extended to transform Matsubara Green functions 



in order to solve the self-consistent equations for the 
BCS-BEC crossover [3(j. Here, three-dimensional spa- 
tial Fourier transformations of isotropic functions and 
discrete Fourier sums with Matsubara frequancies were 
considered. These Fourier transformations are used also 
in the present paper for the numerical calculations. Only 
a few modifications and optimizations have been made 
over the years. The basic principles of the special numer- 
ical Fourier transformation are descibed in the appendix 
of Ref. [3fj| . Here we present the fundamental formulas in 
order to make the numerical method available for appli- 
cations. 

In order to perform a discrete Fourier transformation 
the following sum must be evaluated 



f(k)= J2 Aze^/Or) 



(Bl) 



where a; is a discrete variable with constant step width 
Ax. In this formula and in the following formulas the 
sum over x is defined as a trapezoid sum. This means 
that the first term and the last term in the sum arc multi- 
plied by a factor i, respectively. The continuous Fourier 
transformation is defined by a related integral which is 
obtained from (|B1[) in the limit Ax — > 0. 

We assume that the function values are known in a fi- 
nite subset of points Xj according to f{xj) = aj where 
j = 0, 1, . . . , N. The points Xj cover the whole interval 
between o; m i n and x max on a logarithmic scale so that 
imin = x Q < xi < . . . < x N -i < x N = x max . Conse- 
quently, the Fourier sum (|B1[) can be divided into a sum 
of ./V trapezoid sums according to 



N-l Xj+i 
j=0 x=Xj 



(B2) 



Now, we assume that the function is given by the cubic 
spline polynomial 

f(x) = aj + bj(x - xj) + Cj(x - Xj) 2 + dj(x - x 3 f (B3) 

if x is located in the interval Xj < x < Xj + \. The spline 
coefficients aj, bj, Cj, and dj are calculated numerically. 
Inserting the cubic spline polynomial (|B3[) into the for- 
mula (|B2[) we find that the trapezoid sums within the 
curved brackets can be evaluated exactly. Thus, as a 
result we obtain the Fourier transform 



N-l 



f(k) = ^{ajlf\k) + bjl^(k)+cjl^(k)- 



(3) 



djlj 



3=0 



where 



(*)} 
(B4) 

i)_l]" 
(B5)~ 

By construction a cubic spline function and its first 
two derivatives are continuous. These facts imply the 
following continuity conditions 



&\k) 



ik.i 





a 


-Ax , 








cotl 


'fcAxj 






L 2i 1 
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f{x j+ i) = aj + bj(x j+1 - Xj) + Cj(x j+1 - Xjf + d 3 {x j+1 - Xj f = a j+1 , (B6) 
f'(x j+1 ) = bj + 2cj{x j+1 -xj) + 3dj(x j+1 - Xjf = b j+1 , (B7) 
f"( x j+i) = 2c j + 6d J (x J+1 -Xj) = 2c j+1 , (B8) 



which may be used to regroup the terms in (|B4|) . Conse- 
quently, as a result we obtain the alternative formula 

f(k) = J^(k)[e lkXN a N - e tkx °a ] 

+ J^(k)[e lkXN b N -e lkx °b ] 

+ J^(k)[e lkXN c N -e lkxo c ] (B9) 

N-l 

+ J (3) (fc) [(e lfcEj+1 - e tkx i)dj] 

3=0 

where 

^w- Hint «*(*£)]■ <™> 

The terms with spline coefBcicnts aj , bj , and Cj have can- 
celled for j = 1, 2, . . . , N— 1. In the limit k —» the func- 
tions (|BlU|) diverge according to j(")(fc) - |fc|-( n+1 ). 
For this reason, the alternative formula (|B9[) can be ap- 
plied numerically only for large k so that | kxj | > 1 for all 
j = 0, 1, . . . , N. On the other hand, the functions (|B5|) 
are finite in the limit k — > so that the formula (|B4j) 
can be applied numerically for small fc where | kxj | < 1 
for all j = 0, 1, . . . , N. In practice we use a combination 
of both formulas (|B4|) and (|B9j) . Which formula is used 
for a particular j we decide by considering the value of 
| kxj | and comparing this value with 1 . In this way we 
obtain a special numerical Fourier transformation which 
is stable and reliable for points Xj and ki distributed on 
a logarithmic scale over many decades. 

We have derived our special numerical Fourier trans- 
formation for discrete variables x with a finite constant 
step width Ax. The continuous Fourier transformation is 
obtained easily and naturally by taking the limit Ax — > 
in the functions (|B5[) and (|B10[) which is well defined. 



In order to transform the Green and vertex functions 
forward and backward, we need two kinds of Fourier 
transformations. First we transform between the Mat- 
subara frequencies and the imaginary time variable. In 
this case we can apply a continuous (forward) and a dis- 
crete (backward) Fourier transformation (|F31[) directly. 
Secondly we transfrom between the wave vector and the 
spatial coordinate in d = 3 dimensions. Since the func- 
tions are spherically symmetric, an integration over the 
angles can be performed, so that the resulting transfor- 
mation integrals are one dimensional depending only on 
radial variables, a radial wave number and a radial space 
coordinate, respectively. For d = 3 the transformation 
integrals can be recast into a one-dimensional continu- 
ous Fourier transformation so that our special numerical 
Fourier transformation (|Blj) can be used once again. 

In practice we use N = 300 points for all variables. The 
values of the wave numbers and the values of the space 
coordinates are distributed on logarithmic scales over six 
decades, respectively. The Matsubara frequencies are dis- 
tributed on a logarithmic scale over about twelve decades. 
The imaginary time variables are distributed appropri- 
ately over a finite interval with two logarithmic scales, 
one for each boundary. 

The Green and vertex functions are singular and have 
slowly decaying long tails. For this reason, reference func- 
tions must be subtracted which remove the singularities 
and the long tails. The reference functions are derived 
from free Green functions and the two-particle scattering 
amplitude (T matrix) . For these reference functions ana- 
lytical expressions must be available in all Fourier repre- 
sentations. The difference functions f(x) which are even- 
tually transformed by our numerical method (|B1[) must 
be smooth in x and decay accoording to f(x) ~ x~ 2 or 
faster for \x\ — > oo. 
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